Human population dynamics in Upper Paleolithic Europe inferred from fossil dental phenotypes

Despite extensive archaeological research, our knowledge of the human population history of Upper Paleolithic Europe remains limited, primarily due to the scarce availability and poor molecular preservation of fossil remains. As teeth dominate the fossil record and preserve genetic signatures in their morphology, we compiled a large dataset of 450 dentitions dating between ~47 and 7 thousand years ago (ka), outnumbering existing skeletal and paleogenetic datasets. We tested a range of competing demographic scenarios using a coalescent-based machine learning Approximate Bayesian Computation (ABC) framework that we modified for use with phenotypic data. Mostly in agreement with but also challenging some of the hitherto available evidence, we identified a population turnover in western Europe at ~28 ka, isolates in western and eastern refugia between ~28 and 14.7 ka, and bottlenecks during the Last Glacial Maximum. Methodologically, this study marks the pioneering application of ABC to skeletal phenotypes, paving the way for exciting future research avenues.

Dots indicate instances where pairwise correlations could not be calculated due to insufficient data or variation.Significant associations are denoted with a black star.Significance was evaluated both without (lower triangle) and with (upper triangle) the application of the Benjamini-Hochberg correction for multiple testing.For the meaning of trait abbreviations, refer to Table S1.Out of 1,781 estimated pairwise correlations, 72 remain significant after the Benjamini-Hochberg correction, resulting in an overall ~4% significant trait associations in the dataset.S1.Out of 336 estimated pairwise correlations, two remain significant after the Benjamini-Hochberg correction, resulting in an overall ~0.6% significant trait associations in the dataset.

Fig. S9. Correlations among dental trait expressions in the 'modified key tooth' dataset.
Presented are pairwise tetrachoric correlations, with color coding indicating the direction and strength of the correlation.Dots indicate instances where pairwise correlations could not be calculated due to insufficient data or variation.Significant associations are denoted with a black star.Significance was evaluated both without (lower triangle) and with (upper triangle) the application of the Benjamini-Hochberg correction for multiple testing.For the meaning of trait abbreviations, refer to Table S1.Out of 406 estimated pairwise correlations, six remain significant after the Benjamini-Hochberg correction, resulting in an overall ~1.5% significant trait associations in the dataset.S9.Dental trait sample sizes for the 'key tooth' dataset.Shown are 20 ASUDAS dental non-metric traits, the observable trait sample size (n), the number of traits scored as absent (-), the number of traits scored as present (+), and the frequencies of traits scored as present (%) for six spatiotemporal groups defined for demographic modeling.For the meaning of trait abbreviations, refer to Table S1.'NA' indicates the inability to calculate trait frequencies because the trait could not be observed in the spatiotemporal group.The six groups are defined as follows: West=extending from present-day Portugal to Germany; East=extending from present-day Italy to Western Russia; MPG=Middle Pleniglacial (~47-28 kya); LPG=Late Pleniglacial (~28-14.7 kya); LG&EH=Late Glacial to Early Holocene (~14.7-7kya).

Fig. S2 .
Fig. S2.LDA plot of the 'Discontinuity in the West' set of models using the 'modified key tooth' dataset.Colors correspond to models.The location of the observed data is indicated by a black star.

Fig. S3 .
Fig. S3.LDA plot of the 'Discontinuity in the East' set of models using the 'modified key tooth' dataset.Colors correspond to models.The location of the observed data is indicated by a black star.

Fig. S4 .
Fig. S4.LDA plot of the best-fitting models identified within the 'Continuity', 'Discontinuity in the West', and 'Discontinuity in the East' set of models using the 'modified key tooth' dataset.Colors correspond to models.The location of the observed data is indicated by a black star.

Fig. S5 .
Fig. S5.LDA plot of the simultaneous model comparison using the 'modified key tooth' dataset.Colors correspond to models.The location of the observed data is indicated by a black star.

Fig. S6 .
Fig. S6.LDA plot of the best-fitting models identified with the step-wise and simultaneous model selection approach using the 'modified key tooth' dataset.Colors correspond to models.The location of the observed data is indicated by a black line.

Fig. S7 .
Fig. S7.Correlations among dental trait expressions in the full dataset.Presented are pairwise tetrachoric correlations, with color coding indicating the direction and strength of the correlation.Dots indicate instances where pairwise correlations could not be calculated due to insufficient data or variation.Significant associations are denoted with a black star.Significance was evaluated both without (lower triangle) and with (upper triangle) the application of the Benjamini-Hochberg correction for multiple testing.For the meaning of trait abbreviations, refer to TableS1.Out of 1,781 estimated pairwise correlations, 72 remain significant after the Benjamini-Hochberg correction, resulting in an overall ~4% significant trait associations in the dataset.

Fig. S8 .
Fig. S8.Correlations among dental trait expressions in the 'key tooth' dataset.Presented are pairwise tetrachoric correlations, with color coding indicating the direction and strength of the correlation.Dots indicate instances where pairwise correlations could not be calculated due to insufficient data or variation.Significant associations are denoted with a black star.Significance was evaluated both without (lower triangle) and with (upper triangle) the application of the Benjamini-Hochberg correction for multiple testing.For the meaning of trait abbreviations, refer to TableS1.Out of 336 estimated pairwise correlations, two remain significant after the Benjamini-Hochberg correction, resulting in an overall ~0.6% significant trait associations in the dataset.

Table S1 . Description of dental traits used in this study. Dental
ancestral and presence is derived.'No data' indicates the inability to calculate Pearson correlation due to the absence of published worldwide modern population trait frequencies.

Table S2 . Dental trait sample sizes for the 'modified key tooth' dataset.
Shown are 20 ASUDAS dental non-metric traits, the observable trait sample size (n), the number of traits scored as absent (-), the number of traits scored as present (+), and the frequencies of traits scored as present (%) for six spatiotemporal groups defined for demographic modeling.For the meaning of trait abbreviations, refer to TableS1.'NA' indicates the inability to calculate trait frequencies because the trait could not be observed in the spatiotemporal group.The six groups are defined as follows: West=extending from present-day Portugal to Germany; East=extending from presentday Italy to Western Russia; MPG=Middle Pleniglacial (~47-28 kya); LPG=Late Pleniglacial (~28-14.7 kya); LG&EH=Late Glacial to Early Holocene (~14.7-7kya).

Table S3 . Confusion matrix and model selection results for the 'Continuity' set of models using the 'modified key tooth' dataset.
The most supported model is highlighted in bold.

Table S4 . Confusion matrix and model selection results for the 'Discontinuity in the West' set of models using the 'modified key tooth' dataset.
The most supported model is highlighted in bold.

Table S5 . Confusion matrix and model selection results for the 'Discontinuity in the East' set of models using the 'modified key tooth' dataset.
The most supported model is highlighted in bold.

Table S6 . Confusion matrix and model selection results for the final comparison among the three best-fitting models using the 'modified key tooth' dataset.
The most supported model is highlighted in bold.

Table S7 . Confusion matrix and model selection results for the simultaneous model comparison using the 'modified key tooth' dataset.
The most supported model is highlighted in bold.

Table S8 . Confusion matrix and model selection results for the comparison between the two best-fitting models identified with the step-wise and simultaneous models selection approach using the 'modified key tooth' dataset.
The most supported model is highlighted in bold.